clear all;
clc;
disp('To calculate default prob and realized credit cycle index Z');
load TP06_Y_RJ_fin_R5_G50_250_MR7.mat
% freqj=4;% 4 quarters
year_list=[1997:2009]';
year_list_text=num2str(year_list);
% quarter_list=[1:4]';
% quarter_list_text=num2str(quarter_list);
Rate_list=[0,1:0.5:5];
% default_prob=zeros(length(year_list)*freqj,3);
default_prob=zeros(length(year_list),2);
InvCDF=zeros(length(year_list),2);
Z=zeros(length(year_list),2);
filename='PD et Z_Yr_G50_250_MR7.xls';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate default probability  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:13 %start from 1997 to 2009
%     for j=1:freqj
        default_prob(i,1)=year_list(i,1);
%         default_prob(4*(i-1)+j,1)=year_list(i,1);
%         default_prob(4*(i-1)+j,2)=quarter_list(j,1);
        default_prob(i,2)=(sum(final_result(i,2:10,11))/sum(final_result(i,2:10,13)))/(1-(sum(final_result(i,2:10,12))/sum(final_result(i,2:10,13))));
           %We exclude the debtors whose initial rating is zero           
%      end
end
xlswrite(filename, default_prob, 'default_prob');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To calculate real credit cycle index Z  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:13 %start from 1997 to 2009
%     for j=1:freqj
        InvCDF(i,1)=year_list(i,:);
%         InvCDF(4*(i-1)+j,1)=year_list(i,:);
%         InvCDF(4*(i-1)+j,2)=quarter_list(j,:);
        InvCDF(i,2)=norminv(default_prob(i,2),0,1);
%     end
end
xlswrite(filename, InvCDF, 'InvCDF');

for i=2:13 %start from 1998 to 2009
%     for j=1:freqj
        Z(i,1)=year_list(i,:);
%         Z(4*(i-2)+j,2)=quarter_list(j,:);
%         Zt=InvCDF(4*(i-2)+j,3);
%         miu=mean(InvCDF(5:4*(i-2)+j,3));
%         sig=std(InvCDF(5:4*(i-2)+j,3));
        Zt=InvCDF(i,2);
        miu=mean(InvCDF(2:i,2)); % questions??????? time point
        sig=std(InvCDF(2:i,2));
        Z(i,2)=(-(Zt-miu))/sig; % std is devided by N-1, unbiased           
%      end
end
xlswrite(filename, Z, 'Z');
disp('fin');

% Z is calculated by mean of t time point. Forecast Z is calculated by mean
% of t-1 time point.